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Abstract. When integrating functions that have poles outside the interval of 
integration, but are regular otherwise, it is suggested that the quadrature rule 
in question ought to integrate exactly not only polynomials (if any), but also 
suitable rational functions. The latter are to be chosen so as to match the most 
important poles of the integrand. We describe two methods for generating such 
quadrature rules numerically and report on computational experience with them. 



Introduction 

Traditionally, Gauss quadrature rules are designed to integrate exactly polynomi- 
als of maximum possible degree. This is meaningful for integrand functions that are 
"polynomial-like" . For integrands having poles (outside the interval of integration) it 
would be more natural to include also rational functions among the functions to be 
exactly integrated. In this paper we consider n-point quadrature rules that exactly 
integrate m rational functions (with prescribed location and multiplicity of the poles) 
as well as polynomials of degree 2n — m — 1, where < m < In. The limit case m = 2n, 
in which only rational functions are being integrated exactly, is a rational counterpart 
of the classical Gauss formula; the latter corresponds to the other limit case m = 0. 

In §1 we characterize these new quadrature rules in terms of classical (polynomial) 
Gauss formulae with modified weight functions. We also identify special choices of 
poles that are of interest in applications. The computation of the quadrature rules is 
discussed in §2, and numerical examples are given in §3. 
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1. Gauss quadrature for rational functions 

Let dX be a measure on the real line having finite moments of all orders. Let e C, 
\x = 1, 2, . . . , M, be distinct real or complex numbers such that 



C M ^0 and l + C M i^0 for t G supp(rfA), fi = 1, 2, . . . , M . (1.1) 

For given integers m, n with 1 < m < 2n, we wish to find an n-point quadrature rule 
that integrates exactly (against the measure dX) polynomials of degree 2n — m — 1 as 
well as the m rational functions 

(l + C M t)~ s , /i = l,2,...,M, s = l,2,...,s„, (1.2) 

where > 1 and 

M 

In the extreme case m = 2n (where polynomials of degree -1 are understood to be 
identically zero) the formula integrates exactly 2n rational functions (with poles of 
multiplicities s M at — 1/C/*), but no nontrivial polynomials. The formula, therefore, 
can be thought of as the rational analogue of the classical Gauss formula; the latter 
corresponds to the other limit case m = M = 0. 

The solution of our problem is given by the following theorem. 



THEOREM 1.1. Define 



M 

Wm(*)= nU + CM*)**, (1-4) 



a polynomial of degree m. Assume that the measure dX/uo m admits a (polynomial) 
n-point Gaussian quadrature formula 

L M ^ = £ W »fffi + R "W> R n(*2n-l) = 0, (1.5) 
Jti, UJ m {l) v= i 

with nodes t^j contained in the support of dX, 

t„ e supp(rfA). (1.6) 

Define 

t v = t G v , \ v = w°w m (t% v = 1,2,..., n. (1.7) 

Then 

„ n 

g(t)dX(t) = ^ X v g(U) + R n (g), (1.8) 



where 



g(t) = (l + (^t)- s , fjL = l,2,...,M; s = 1, 2, . . . , s M , 
/?„(#) = /./' { (1.9) 

g G P2n-m-l • 



Conversely, (1.8) wrai/i ^ G supp(dA) and (1.9) imply (1.5), (1.6) with t^,w^ as 
defined in (1.7). 

Remark. Theorem 1.1, for real and either all s^ = 1 and to = 2n, or all but one 
s^ = 2 and m = 2ra — 1, is due to Van Assche and Vanherwegen [13]. The quadrature 
rule (1.5), especially its convergence properties for analytic functions /, has previously 
been studied by Lopez and Illan [9, 10]. 

Proof of Theorem 1.1. Assume first (1.5), (1.6). For p, — 1, 2, . . . , M; s — 1, 2, . . . , s M , 
define 

q (f) = ^ . (1.10) 
Since m < 2n and s > 1, we have g MjS G P m -s C P2n-i, and therefore, by (1.5), 

" (i + o^) s ~ (i + ^) s ' 

where (1.7) has been used in the last step and none of the denominators on the far 
right vanishes by (1.6) and (1.1). This proves the assertion in the top line of (1.9). The 
bottom part of (1.9) follows similarly: Let p be an arbitrary polynomial in P2 n - m -i- 
Then, since poJ m G P2n-i, again by (1.5) and (1.7), 

p(t)d\(t) = [p(t)u; m (t) dMn 



R JR uJ m (t) 

n n 
v=l v=\ 



To prove the converse, we first note that is well defined by (1.7), since u) m (t v ) ^ 
by the assumption on t v and (1.1). One then easily verifies that (1.5) holds for all 
polynomials (1.10) (of degree < to) and all polynomials of the form puo m where 
P G P2n-i- m . The collection of these polynomials, however, spans P2n-i- n 



We will concentrate on six special choices of the parameters C M that are of interest 
in applications. 

Case 1 (Simple real poles). All = 1 in (1.2) (hence M = to), and all are real, 
distinct, and nonzero, 

C„ = &eR, 6^0, u = 1,2,..., to. (1.11a) 
In this case the polynomial u m has the form 

m 

Wm(<)= 11(1 + ^6R. (1.11b) 

If the support of d\ is an interval, u m does not change sign on it because of (1.1). 

Case 2 (Simple conjugate complex poles). All s^ — l (hence M = m), m even, and 
the occur in to/2 (distinct) pairs of conjugate complex numbers (cf. [9]), 

Cv = iu + iVu, Cu+ m /2 = iv-'iVv, v = 1,2,..., m/2, (1.12a) 
where £„ G R and ^ > 0. Here, 

m/2 

"m(f)=II[( 1 + &*) 2 + '# 2 ], (l-12b) 

I/=l 

which is strictly positive for real t. 

Case 2' (Simple conjugate complex poles plus a simple real pole). All s^ — l (hence 
M — to), to (odd) > 3, and, slightly changing the indexing of the £'s, 

Co e R, C = 6 + ^, C+(m-i)/2 = £„ ~ iVv, ^ = 1,2, . . . , (m - l)/2, (1.13a) 
where Co = £0 7^ and (,6R, to > for 1 < v < (to - l)/2. Then 

(m-l)/2 

w m (i) = (l + £ f) II [(l + ^) 2 + ^* 2 ]- (l-13b) 
i/=i 



Case 3 (Real poles of order 2). All s^ = 2 in (1.2) (hence 2M = to), and all C^ are 
nonzero, real and distinct, 

C, = &,GR, £,^0, s„ = 2, z/= 1,2,..., to/2. (1.14a) 



The polynomial u m now has the form 

m/2 



LO, 



»(o = n( i +^*) 2 (i-i4b) 



and is nonnegative for real t, and positive on the support of dX. 

Case 3' (Real poles of order 2 plus a simple real pole). Here, all = ^ are nonzero, 
real and distinct, s v — 2 for v — 1, 2, . . . , M — 1 and s M = 1- Thus, m = 2M — 1, and 



M - x m + 1 



in 



(t) = (1 + £ M t) n (1 + M = — ^— , m(odd) > 3. (1.15) 



i/=i 



If M = n, i.e., m = 2n — 1, the quadrature rule (1.5) is then identical with the 
"orthogonal quadrature rule" of [13], having as nodes the zeros of the rational function 
(1 + C,nt)~ l + Z)"=i Qi/(1 + C^)^ 1 which is orthogonal (relative to the measure dX) to 1 
and to (1 + Cm^) _1 ' A* = 1, 2, . . . , n — 1. As in Case 1, the polynomial u m preserves its 
sign on the interval on which d\ is supported. 

Case 4 (Conjugate complex poles of order 2). All s M = 2 (hence 2M = m), m = 
(mod 4), and the fj, — 1, 2, . . . , m/2, occur in m/4 conjugate complex pairs, similarly 
as in Case 2. Thus, 

m/4 

"m(*)= Il[(l + ^) 2 + ^ 2 ] 2 - (1-16) 

l/=l 

In all six cases, the measure d\/u m admits a Gaussian n-point formula for each 
n — 1, 2, 3, . . . , so that the assumption of Theorem 1.1 is fulfilled for each n. 
Putting f(t) = uj m {t)g{t) in (1.5) and using (1.7), we get 

„ n 

g(t)d\(t) = £ X u g(t u ) + R°(u m g), (1.17) 

where from the well-known expression for the remainder term of Gaussian quadrature 
rules, one has 

nSM = 7nM) (2n) (r), in = ^ 2n) f n ■ (i-i8) 

Here, r is some number in the smallest interval containing the support of dX, and 
Pk — (3k(dX/uj m ) are the /5-recursion coefficients for the measure dX/uj m (cf. (2.1) 
below). The latter are computed as part of the algorithms to be described in the next 
section. 



2. Computation of the quadrature rule (1.5) 



We propose essentially two methods for generating the basic quadrature rule (1.5), 
the first being most appropriate if the support of dX is a finite interval, the other more 
effective, though possibly slower, when the support interval of dX is unbounded. 

2.1. Method based on partial fraction decomposition and modification 
algorithms. To compute the n-point formula (1.5), it suffices to compute the nth- 
degree orthogonal polynomial if n ( ■ ) = ir n ( ■ ; dX) relative to the measure dX = dX/uj m , 
or, more precisely, the recursion coefficients a k = a k (dX), j3 k — Pk(dX), k = 0, 1, ... , 
n — 1, in the three-term recurrence relation satisfied by these (monic) polynomials: 



fc = 0,l,...,n-l, (2- 1 ) 



7T (t) = l, 7r_!(t) = 0. 

The nodes t„ and weights w„ in (1.5) can then be obtained by standard techniques 
via an eigensystem problem for the (symmetric, tridiagonal) Jacobi matrix of order n 

having the a k , k — 0, 1, . . . , n — 1, on the diagonal, and \f$k, k — 1, 2, . . . , n — 1, on 
the side diagonals (see, e.g., [7], [4, §6]). The coefficients a k , f3 k in turn are expressible 
in terms of the orthogonal polynomials n k as 

(tn k ,n k ) 

ak = -r„ — r^r , < k < n - 1, 

{^k, ^k) ^.2) 

A) = (tto, tto), $k = ,S nk,n ^ \ , 1 < k < n - 1, 
where ( • , • ) denotes the inner product 



(u,v) = / u(t)u(t)d\(t). (2.3) 

(If the error constant 7„ in (1.18) is desired, one needs to compute, in addition, f3 n .) 

The basic idea of computing the coefficients in (2.2) is as follows. Suppose we can 
construct an iV-point quadrature rule for dX = dX/u m , where N > n, which is exact 
for polynomials of degree < 2n — 1: 



p(t)dX(t) = jrw kP (T k ), peP 2n _!. (2.4) 
R k=i 

Here the weights W k are not necessarily all positive. Denote the discrete measure 
implied by the sum on the right by dA N : 

[ P (t)dA N (t) = jrw kP (T k ). (2.5) 
JR k=i 



Then from the formulae in (2.2) one easily sees by induction that 



a k (dX) = a k (dA N ), 

fc = 0,l,...,n-l. (2.6) 

p k (d\) = p k (dA N ), 

Thus, the desired recursion coefficients are the first n of the a- and /3-coefficients 
belonging to the discrete measure dA N . These can be generated by Stieltjes's procedure, 
which is implemented in the routine sti of [4]. (The faster routine lancz of [4], 
implementing the Lanczos method, would also be applicable here, even though dA N is 
not necessarily a positive measure.) 

We next show how a quadrature rule of type (2.4), with N = 0(mn), can be 
constructed by means of partial fraction decomposition and suitable modification al- 
gorithms. For this, we consider separately Cases 1-3' identified in §1. The analysis of 
Case 4 becomes so tedious that we will not pursue it any further in this context; see, 
however, §2.2. 

2.1.1. Simple real poles. We set up the partial fraction decomposition of l/u m in 
the form 



1 1 



E TT~77TF~\ > ( 2 - 7 ) 



where 



and an empty product in (2.8) (when m — 1) is to be taken as 1. Then, with dX ■ 

dX/ u m , 



The integrals on the right involve measures c u dX modified by linear divisors. For such 
measures, the associated recursion coefficients can be obtained from those of c v dX 
(assumed known) by a suitable modification algorithm (cf. [4, §5]). Unless x v = — 1/£„ 
is very close to the support interval of dX, the most appropriate algorithm is the 
one embodied in the routine gchri of [4] with iopt = 1. Otherwise, the routine 
chri of [4] (again with iopt = 1) is preferable. A basic ingredient of the routine 
gchri is the modified Chebyshev algorithm (cf. [2, §2.4]) using modified moments 
Jr 7Tfc(t; dX)c v dX(t) / (t — x u ), k — 0, 1, 2, . . . , 2n — 1. These in turn are generated by 
backward recurrence as minimal solution of the three-term recurrence relation for the 
measure dX; cf. [1, §5]. 



Having obtained, in whichever way, the first n of the a- and /^-coefficients for the 
modified measure c u dX(t)/(t — x u ), and hence the Gaussian quadrature formula^ 



n 



■ K P(t) = E/rM4% V e P 2n -u (2.9) 

via eigensystem techniques, we then get 

d\(t) ™ f c v dX(t) 



m n 



u=l r=l 

hence the desired quadrature rule (2.4), with N = mn and 

Ti i\ , = 

± (y— l)n+r v r j 

u= 1,2, ...,m; r = l,2,...,n. (2.10) 

W(„_l)n+r = lt;M, 

The procedure described works best if the support of g?A is a finite interval. Other- 
wise, the modified Chebyshev algorithm underlying the procedure is likely to suffer 
from ill-conditioning; cf. Example 3.4. Another difficulty that may adversely affect the 
accuracy of the results, in particular if m = 2n, is the possibility that the constants 
c u sgn tgsU pp( dA ) (t + become very large and alternate in sign; cf. Example 3.2. 

This will cause serious cancellation errors in evaluating inner products relative to the 
measure dAjy (there being blocks of weights Wk which are very large positive alternating 
with blocks of weights which are very large negative). In such cases, either m has to be 
lowered, perhaps down torn = 1, or else the method discussed in §2.2 invoked, which 
will be more effective (but possibly more expensive). 

2.1.2. Simple conjugate complex poles. We now consider Case 2 of §1, i.e., conjugate 
complex parameters = £„ + ir] u , ( u+m /2 = (v, where £ v G R, rj v > and m is even. 
In this case, an elementary computation yields the partial fraction decomposition 



m/2 



+ d v t 



t 6 R, (2.11) 



^n order to produce positive /^-coefficients, as required in the routine for Gauss quadrature for- 
mulae, one inputs the measure \c u /(t — x v )\dX(t) and, if this entails a change of sign, reverses the sign 
of all Gauss weights after exiting from the Gauss quadrature routine. 



where 



c - = i Im ^ + sfe Re ^) ' (2 12) 

and 

m/2 + in ) 2 

with pi — 1 if m — 2. One can then proceed as in §2.1.1, except that the modification 
of the measure d\ now involves multiplication by a nonconstant linear function (if 
d v 7^ 0) in addition to division by a quadratic. The former modification is handled 
by the routine chri of [4] with iopt = 1, the latter by the routine gchri with iopt 
= 2 (or, if more appropriate, by chri with iopt = 5). The quadrature rule (2.4) so 
obtained has N = mn/2. 

If the poles — 1/C/u are located in conjugate pairs on a line parallel to the imaginary 
axis, then by an elementary calculation one can show that all p v are real, hence d v = 0, 
and there is no need to call chri. 

2.1.2'. Simple conjugate complex poles plus a simple real pole. We are now in Case 
2' of §1, with m odd, Co = £o £ R- an d the remaining conjugate complex as in Case 
2. This yields 

-^tt = rV^ + ( "^ iT C 1±^A t £ R (2.14) 

where 



c, 



/■m— 2 
SO 



ni,=7 1)/2 [(6-o 2 +^ 



Im P :+ — ^RepM, (2.15) 



d' v = — Im p' v 



and 



4v - so + ^ 

with ]9j, the same as in (2.13) with m replaced by m — 1. The technique called for is 
a combination of the one in §2.1.1, to deal with the first term in (2.14), and the one 
in §2.1.2, to deal with the remaining terms, and yields a quadrature rule (2.4) with 
N = (m + l)ra/2. 



2.1.3. Real poles of order 2. This is Case 3 of §1, and leads to the partial fraction 
decomposition 

j m/2 / d \ 

M*) = S v^TT/^ + (t + i/o 2 J ' (2 ' 17) 

otm-3 g M 

c „ = —J±i (2.18) 

cm— 4 

d u = , (2.19) 

where c± — 0, d± — £f 2 when m — 2. Here, A" = mn in (2.4). 

2.1.3'. Real poles of order 2 plus a simple real pole. Similarly as in §2.1.3, the partial 
fraction decomposition has now the form 

1 r' M ~ l ( r> rl' \ 



(2.20) 



c 



fm-2 

/ _ >M 
M — 



km + 2(^- Cm) E^T 1 t%- 



cm— 4 



& - Cm) ntLT 1 ^ - 4>) 5 



Empty sums and products (when M — 2) have their conventional values and 1, 
respectively. Again, A" = mn in (2.4). 

The presence of two terms in the summations of (2.17) and (2.20) complicates mat- 
ters considerably, as they call for two applications of the routine gchri: First, we must 
generate sufficiently many of the recursion coefficients for the measure d\(t)/(t — x v ), 
x v = — l/£i/, in order next to generate the desired recursion coefficients for d\(t)/(t — 
x u ) 2 by backward recursion - a recursion based on the recurrence relation generated 
in the first application of gchri (which in turn requires backward recursion!). The 
procedure nevertheless works well if the x v are not too close to the support interval of 
dX; see Example 3.3. 



2.2. Discretization method. In this method, the inner product (2.3) is approx- 
imated by a discrete (positive) inner product, 

(u, v) = [ u(t)v(t) ^§ « £ c^rf } ) =: («, JV > n, (2.21) 



whereupon the formulae (2.2) are applied with the inner product ( • , • ) replaced by 
( • , • )jv throughout. This yields approximations 

«fc,w~«fc, Pk,N~Pk, fc = 0,l,...,ra-l. (2.22) 

In effect we are generating the polynomials orthogonal with respect to the discrete 
inner product ( • , • ) jv in order to approximate the desired orthogonal polynomials. 

The computation of the approximate coefficients (2.22) can be done by either 
Stieltjes's procedure or Lanczos's algorithm (cf., e.g., [3, §§6-7]). Both are imple- 
mented in the routine mcdis of [4]. 

With any reasonable choice of the discretization (2.21), it will be true that the 
procedure converges as iV — > oo, 

lim a KN = d k , lim f3 KN = $ k , < k < n - 1. (2.23) 



A natural choice, indeed, is given by 

rf } = A N \d\), u[ N) = - /,• =1.2 V. (2.2 I) 



w i N \d\) 
t 



where t k N \d\) are the zeros of the orthogonal polynomial ttn{ ■ ]d\), and w[ N \d\) 
the respective Christoffel numbers. 

The discretization method is conceptually simpler, and sometimes more stable, than 
the methods of §2.1, but may become significantly more expensive, regardless of the 
choice of to, if poles are close to the interval of integration, or if high accuracy is 
desired; cf. Examples 3.1 and 3.5. Note also that Case 4 that was skipped in §2.1 can 
easily be handled by the present method; see Example 3.6. 

3. Numerical Examples 

All examples in this section were computed on the Cyber 205 in both single and 
double precision. The respective machine precisions are 7.11 x 1CT 15 and 5.05 x 10~ 29 . 



Example 3.1. h{u>) = j\ dt, u > 1. 

Here, d\(t) = dt, and the poles of the integrand are located at the integer mul- 
tiples of uj. It is natural, then, to make our quadrature rule (1.8) exact for to ele- 
mentary rational functions matching the to poles closest to the origin, say those at 



— (m/2)u!, . . . , — uj, uj, . . . , {m/2)uj when m is even. This suggests to identify — 1/C/i hi 
(1.2) with these poles, i.e., in (1.11a) to set 

£, = (-l)7(c4(z/ + l)/2j), v = l,2,...,m. (3.1) 

Best accuracy is expected when m = 2n, in which case the method described in §2.1.1 
was found to work rather well, the only difficulty being the relatively slow conver- 
gence of the backward recurrence algorithm for computing the 2n modified (Legendre) 
moments associated with the measure dt/(t ±uj) when uj is very close to 1. For single- 
precision accuracy e = \ x 1CT 10 and double-precision accuracy e d = \ x 1CT 25 , and 
n = 20, the respective starting indices ko and k$ in the backward recursion yielding 
the desired accuracy are shown in Table 3.1 for selected values of uj. 





h 


Ud 


2.0 


50 


63 


1.5 


53 


71 


1.1 


67 


106 


1.01 


124 


247 



TABLE 3.1. Starting indices for backward 
recurrence when n = 20 

Other than that, the method appears to be very stable and produces quadrature 
rules that are rapidly converging. In Table 3.2, the results of the n-point rule (1.17) 



UJ 


n 


n-point rational Gauss 


In 




err. 


Gauss 


2.0 


1 


2.1 


3.94(- 


-1) 


1.43 






4 


2.33248722 


3.50(- 


-7) 


7.18 


:-s) 




7 


2.332487232246550235 


2.61(- 


-15) 


2.73 


:-s) 




10 


2.332487232246550241107076 


1.48(- 


-24) 


1.02 


:-n) 


1.1 


2 


4.43 


1.73(- 


-2) 


2.60 






5 


4.467773637 


2.00(- 


-9) 


2.09 


:-2) 




8 


4.46777364638776571 


5.61(- 


-18) 


1.53 


:-3) 




11 


4.467773646387765789236123 


1.66(- 


-27) 


1.09 


:-4) 


1.01 


3 


8.429 


2.53(- 


-4) 


4.20 






6 


8.4301845803 


6.27(- 


-12) 


1.85 






9 


8.4301845804708420582 


7.52(- 


-21) 


8.37 


:-2) 




12 


8.430184580470842058971264 


1.23(- 


-30) 


3.75 


:-2) 



TABLE 3.2. Numerical results for I\{uS), error constants, 
and comparison with Gauss quadrature 

applied to g(t) = {nt/uj)/ sm.^t/uj) in double precision are shown for uj = 2, 1.1 and 
1.01, along with the error constants 7„ of (1.18). Also shown in the last column are 
the relative errors of the n-point Gauss-Legendre rule. For uj — 2, the exact answer is 
known to be 8C/7T, where C is Catalan's constant (cf. [8, Eq. 3.747(2)]). The value 
shown in Table 3.2 for n — 10 agrees with it to all 25 decimal digits given. Ordinary 
Gauss-Legendre quadrature is seen to converge rather slowly, as uj approaches 1. In 
contrast, convergence of the rational Gauss quadrature rule is fast even for uj very close 
to 1. The extra effort required in this case is expended, as illustrated in Table 3.1, at 
the time when the rule is generated. 

method of §2.1 method of §2.2 



UJ 


n 


SP 


DP 


SP 


DP 


2.0 


1 


.001 


.004 


.007 


.178 




4 


.008 


.036 


.010 


.224 




7 


.027 


.125 


.028 


.220 




10 


.065 


.296 


.016 


.423 


1.1 


2 


.002 


.014 


.060 


1.554 




5 


.013 


.063 


.068 


1.649 




8 


.036 


.177 


.098 


1.274 




11 


.080 


.376 


.104 


1.461 


1.01 


3 


.005 


.037 


.460 


20.10 




6 


.020 


.104 


.446 


41.51 




9 


.050 


.244 


.458 


10.44 




12 


.102 


.487 


.525 


12.65 



TABLE 3.3. Timings (in seconds) of the methods in §§2.1 
and 2.2 applied to Example 3.1 

Identical results were obtained by the discretization method of §2.2, but with sub- 
stantially greater effort, particularly for higher accuracies and for uj close to 1. Respec- 
tive timings are shown in Table 3.3, both for single-precision (SP) and double-precision 
(DP) accuracy requirements of | x 10~ 10 and | x 10 -25 , respectively. 

While the choice m = 2n indeed gives best accuracy, other choices of m may be 
preferable if the effort and time to generate the quadrature rule is of any importance. 
It turns out that with the method of partial fractions, m = 2[(n + l)/2j gives almost 
the same accuracy at about half the effort, whereas m = 2 gives considerably less 
accuracy but requires only about one-tenth the effort. The discretization method of 
§2.2, on the other hand, requires essentially the same effort regardless of the choice of 



to. Some timings required to generate the quadrature rules for various to and n, and 
the relative errors achieved, are shown in Table 3.4 



method of §2.1 method of §2.2 
uj n to SP DP SP DP err 



2.0 10 


2 


.008 


.037 


.014 


.399 


1.10(- 


-17) 




10 


.033 


.153 


.015 


.409 


1.47(- 


-25) 




20 


.065 


.296 


.016 


.423 


1.58(- 


-25) 


1.1 11 


2 


.009 


.047 


.095 


1.402 


2.20(- 


-13) 




12 


.045 


.213 


.100 


1.426 


2.80(- 


-23) 




22 


.080 


.376 


.104 


1.461 


6.68(- 


-26) 


1.01 12 


2 


.011 


.065 


.526 


12.65 


1.15(- 


-13) 




12 


.052 


.258 


.511 


12.46 


9.10(- 


-25) 




24 


.102 


.487 


.526 


12.65 


3.55(- 


-27) 



TABLE 3.4. Timings and errors for selected to < 2n 



If to = 2, the values of n for which full accuracy of about 10" 25 is attained are 15, 
21 and 22 for uj = 2.0, 1.1 and 1.01, respectively. Interestingly, the timings involved 
are only about half those for to = 2[(n + 1)/2J shown in Table 3.4. 

Example 3.2. I 2 (u) = J* dt, < u < 1. 

Here we take d\(t) = t" 1 / 2 dt on [0,1]. If we wish to match the first 2n — 1 poles of 
the gamma function at the negative integers as well as the pole at — uj, we set to = 2n 
in 



= ~> f — 2, 3, ... , to. 

The rational n-point Gauss rule (1.7), (1.17), generated by the method of §2.1.1, then 
produces (in double precision) results as shown in Table 3.5, where uj — \. In the last 
column we list the absolute value of the difference between double-precision and single- 
precision results. In contrast to Example 3.1, we now see a case in which the accuracy 
reaches a limit (at about n — 10) and deteriorates, rather than improves, as n is further 
increased. (When n = 20, the calculation even breaks down in single precision!). The 
last column in Table 3.5 provides a clear hint as to what is happening: a steady growth 
in numerical instability. Closer examination reveals the true cause of this instability. 



The constants c v in the partial fraction decomposition (2.7) become very large and 
alternate in sign. Thus, for example, c 18 = —2.3375 . . . x 10 9 and c 19 = 2.3336 . . . x 10 9 
when n = 18. This produces blocks of large coefficients Wk in (2.10) that alternate 
in sign from block to block, causing severe cancellation errors in summations such as 
(2.5) (which are abundant in Stieltjes's algorithm). The phenomenon evidently is a 
manifestation of the asymmetric distribution of the poles of the gamma function. 



n 


n-point rational Gauss 


\DP - SP 


2 


1.746 


2.24(-13) 


6 


1.75012059121 


7.05(-ll) 


10 


1.750120591261335415386 


3.36(-8) 


14 


1.7501205912613354159 


1.52(-5) 


18 


1.7501205912613356 


4.31 (-3) 



TABLE 3.5. Numerical results for hi^o), cu — .5 



The method of §2.2, in contrast, does not suffer from any numerical instability and 
produces for n = 12, with comparable effort, the value 

J 2 (.5) = 1.750120591261335415394610, (3.3) 

believed to be correct to all 25 digits shown. 

Matching only n poles, and thus taking m = n in (3.2), stabilizes the procedure 
of §2.1.1 considerably, and as a consequence produces the correct result (3.3) (except 
for a discrepancy of 1 unit in the last decimal place) for n = 11. An even more stable 
procedure results from taking m = 2 and, amazingly, yields the correct answer (to all 
digits shown!) already for n = 13. 

Example 3.3. I 3 (u) = j\ ( si ^ M Y dt, u > I. 
Similarly as in Example 3.1, we take 

6, = (-l)7K(^+l)/2j), */ = l,2,...,m/2. (3.4) 

We applied the procedure described in §2.1.3 for uj = 2, 1.5, 1.1 and 1.01, both in single 
and double precision, requesting accuracies of e = \ x 10~ 10 and e d = | x 10~ 25 , 





h 


kf 


k 2 


b d 
h 2 


2.0 


130 


163 


50 


63 


1.5 


133 


191 


53 


71 


1.1 


177 


306 


67 


106 


1.01 


384 


727 


124 


247 



TABLE 3.6. Starting indices for two backward 
recurrences when n = 20 

respectively. When m = 2n and n = 20, starting indices k±, kf in the first application of 
the backward recursion that were found to meet the accuracy requirements for the poles 
closest to [-1,1], and the analogous starting indices &2, k$ in the second application, 
are shown in Table 3.6 for the four values of uj. As expected, the procedure becomes 
laborious as uj approaches 1. Selected double-precision results produced by the n-point 
rational Gauss rule, along with error constants, are shown in Table 3.7. The last 
column shows the relative error of results generated by the n-point Gauss-Legendre 
rule. For uj — 2, the exact answer is known to be J 3 (2) = 4 In 2 ([8, Eq. 3.837(2)]) 
and is correctly reproduced to 25 digits when n — 11. Note again the fast convergence 
of the rational Gauss quadrature rule, even for uj very close to 1, in contrast to the 
relatively slow convergence of the ordinary Gauss rule, especially for uj close to 1. 



UJ 


n 


n-point rational Gauss 


In 




err. 


Gauss 


2.0 


2 


2.75 


9.90(- 


-3) 


4.36 


:-2) 




5 


2.77258868 


1.16(- 


-9) 


4.70 






8 


2.7725887222397811 


3.28(- 


-18) 


2.92 






11 


2.772588722239781237668928 


9.72(- 


-28) 


1.52 


:-n) 


1.1 


2 


15.5 


3.00(- 


-2) 


6.69 






6 


16.5328175 


8.54(- 


-12) 


6.79 


:-2) 




10 


16.5328177384604181 


7.21(- 


-24) 


3.42 


:-3) 




14 


16.53281773846041830155898 


2.35(- 


-37) 


1.40 


:-4) 


1.01 


2 


184. 


6.34(- 


-2) 


9.64 






6 


188.674782 


2.20(- 


-11) 


7.06 






10 


188.674784224994172 


1.88(- 


-23) 


4.00 






14 


188.6747842249941742708325 


6.15(- 


-37) 


1.92 





TABLE 3.7. Numerical results for , error constants, and comparison with Gauss 

quadrature 



While the choice m = 2|_(n + 1)/2J produced similar advantages as in Example 3.1 
- an increase of speed by a factor of about 2 at only a slight loss of accuracy — the 
choice m = 2 offered no significant gains in accuracy over the Gauss-Legendre rule, 
unlike m — 4, which did (since a symmetric pair of double poles is now accounted for). 

We also applied the discretization method of §2.2 and obtained identical results 
with somewhat less effort in the case uj — 2, and about the same effort in the case 
uj = 1.1. For uj = 1.01, however, we were unable to attain the requested double- 
precision accuracy with a discretization parameter < 800 (in (2.21)). 



Example 3.4. I 4 = J °° e~ l dt. 

The appropriate measure here is d\(t) = e~ l dt on [0,oo]. Since the integrand has 
poles at the integer multiples of 27ri, we let ( u = —l/(2vm) = i/(2i/n), and thus in 
(1.12) take 

6 = 0, Vu= 7T-, u = l,2,...,m/2. (3.5) 
Zun 

The quantity p v in (2.13) being real, and thus d v = in (2.12), there is no nonconstant 
linear factor in the numerators of (2.11). This simplifies somewhat the procedure in 
§2.1.2, as it obviates the need to apply the routine chri. 

In Table 3.8 we compare the performance (in double precision and for m = 2n) 
of our rational quadrature routine with Gauss-Laguerre quadrature (applied to f(t) = 
t/(e t — 1)) and the Gaussian quadrature rule (applied to f(t) = e~*) associated with 
"Einstein's weight function" t/(e t — 1); for the latter see [6]. The respective relative 
errors are shown in the last two columns. It can be seen that the Gauss-Laguerre and 
Gauss-Einstein quadratures are comparable in accuracy, the former being somewhat 
more accurate for small values of n, the latter for larger values of n. Both quadrature 
rules, however, are incomparably inferior to the rational Gauss formula, which for 
n — 15 produces the true value of the integral, £(2) — 1 = (7r 2 /6) — 1, to 25 correct 
decimal digits. (Actually, the last digit is off by one unit.) The results become even 
slightly more accurate when we choose m = 2|_(n + l)/2j , and are still better, by several 
orders of magnitude, than those for Gauss-Laguerre and Gauss-Einstein quadrature 
when m = 2. 



n n-point rational Gauss err GL err GE 

1 .59 9.76(-2) 4.09(-l) 

5 .644934055 1.50(-5) 2.97(-4) 

10 .644934066848226428 2.22(-8) 1.15(-8) 

15 .6449340668482264364724151 1.59(-11) 3.25(-13) 

TABLE 3.8. Numerical results for I4 and comparison with 
Gauss-Laguerre and Gauss-Einstein quadrature 



The high accuracy of our rational quadrature rules in this example is all the more 
remarkable as the routine gchri, used in their construction (by the methods of §2.1.2), 
is subject to ill-conditioning, causing the recursion coefficients for the relevant orthog- 
onal polynomials to gradually lose accuracy (by as much as 10 decimals, when n — 15 
and m = 2n). 

This weakness is accentuated when one tries to deal with more difficult integrals, 
for example, 

J W = / ~f — 7 Jl + h0tdt= , 9 > 0, (3.6) 

JO 6 — 1 JO 



which has an additional branch point singularity at t — —2/9. Here, when 9 = .75, the 
n-point rational Gauss formula (in double precision and for m = 2n) gives only about 
13 correct decimal places for n = 15, and 18 for n = 30. By the time n reaches 33, 
the ill-conditioning in the routine gchri has built up to such a level that the method 
fails (by producing a negative /^-recursion coefficient). To get higher accuracy, one 
needs to apply the discretization method of §2.2, which is more stable, but becomes 
fairly expensive if pushed much beyond n = 30. Using d\(t) = e _t dt, and hence the 
Appoint Gauss-Laguerre formula, to effect the discretization in (2.21), and requesting 
an accuracy of | x 10~ 25 for the desired recursion coefficients, we have observed timings 
of the order 12-16 seconds, and discretization parameters iV as large as iV = 370, for 
33 < n < 40. The rational Gauss formula so produced then yields relative errors of 
4.26 x 10" 21 for n = 35, and 9.65 x 10~ 23 for n = 40. This is still better, by about 
4 decimal orders of accuracy, than Gauss-Laguerre quadrature applied to the second 
form of the integral in (3.6), and Gauss-Einstein quadrature applied to the first form. 

Generalized Fermi-Dirac integrals (cf. [12]) are similar to 1(9) except that t/(e t — l) 
is replaced by t k /(e~ v+t + 1), where rj is a real parameter and A; = 1/2, 3/2 or 5/2. The 
poles are now located at r\ ± (2u — l)in, v — 1, 2, 3, . . . . The use of rational Gauss 
quadrature to compute such integrals is dealt with elsewhere [5]. 

Example 3.5. I 5 (rj) = / °° —^-^e^dt, r] < 0. 

Again, we take d\(t) = e~ l dt and note that the poles are now at r\ ± 2uni, v = 
0, 1, 2, . . . . Accordingly, in (1.13) we take 

& = --, £v = — 2 , \ 2 2 , Vv = 2 pj~2 2 ' ^= 1,2, ...,(m- l)/2, (3.7) 

and use the procedure of §2.1.2'. Selected results (for m = 2n — 1), comparing rational 
Gauss formulae with Gauss-Laguerre formulae, are shown in Table 3.9. In the case 



V 


n 


n-point rational Gauss 


err 


GL 


- .1 


3 


.4503 


1.16{ 


-1) 




6 


.4501936153 


5.13{ 


-2) 




9 


.450193614441350 


2.70( 


-2) 




12 


.45019361444134784096 


1.55 


:-2) 


-1.0 


2 


.113 


2.14{ 






6 


.1111093520 


5.07{ 


-3) 




11 


.1111093516052317322 


1.81 


:-4) 




16 


.1111093516052317320105065 


1.26 




-10.0 


2 


.122 (-4) 


1.79( 






6 


.113502121(-4) 


1.57( 


-4) 




11 


.1 135021 14635390578(-4) 


7.31 


:-9) 




16 


.1135021146353905701870968(-4) 


1.20 


:-i2) 



TABLE 3.9. Numerical results for 1 5 and comparison with Gauss-Laguerre quadrature 



r] = — .1, we were able to go only up to n = 13; when n = 14, our procedure failed 
by producing a negative /^-coefficient in (2.2). The difficulty is caused by the ill- 
conditioning (mentioned after (2.10)) affecting the modified Chebyshev procedure. 
Even though our procedure was successful for n = 13, it had to work hard to take 
care of the pole at n = — .1: Backward recursion to compute modified moments had 
to start at v — 584 to get single-precision accuracy \ x 1CT 10 , and at v — 2650 to get 
double-precision accuracy | x 10~ 25 . 

Replacing the numerator te~ l in the integrand by t k , where k = 1/2,3/2 or 5/2, 
and adding a factor \Jl + \9t as in (3.6), produces the Bose-Einstein integral whose 
computation by rational Gauss quadrature is discussed in [5]. 

Example 3.6. J 6 = / °° (-^f e^dt. 

Here, as in Example 3.4, we take d\(t) = e~ f dt and parameters £„, i] u as in (3.5), 
except that there are only m/4 of them, m being divisible by 4. In Table 3.10 we give 



n n-point rational Gauss err GL err GE 

2 .47 3.71(-2) 1.61(-2) 

8 .4816405209 1.16(-6) 5.99(-10) 

14 .4816405210580757311 4.36(-9) 7.26(-18) 

20 .4816405210580757313458777 2.80(-ll) 1.09(-25) 

TABLE 3.10. Numerical results for 1$ and comparison with Gauss-Laguerre and 

Gauss- Einstein quadrature 



the results for m = In in even) obtained by the discretization method of §2.2, analogous 
to those of Table 3.8 but using the square of the Einstein function as weight function 
in GE. (The method of §2.1, as mentioned earlier, was not implemented.) What is 
remarkable in this example is the competitiveness of the Gauss-Einstein quadrature 
rule vis-a-vis the rational Gauss rule. 
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